Reactive force field for lithium-aluminum 
silicates with applications to eucryptite phases 
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We have parameterized a reactive force field (ReaxFF) for lithium aluminum silicates using den- 
sity functional theory (DFT) calculations of structural properties of a number of bulk phase oxides, 
silicates, and aluminates, as well as of several representative clusters. The force field parameters 
optimized in this study were found to predict lattice parameters and heats of formation of selected 
condensed phases in excellent agreement with previous DFT calculations and with experiments. We 
have used the newly developed force-field to study the eucryptite phases in terms of their thermo- 
dynamic stability and their elastic properties. We have found that (a) these ReaxFF parameters 
predict the correct order of stability of the three crystalline polymorphs of eucryptite, a, ft, and 
7, and (b) that upon indentation, a new phase appears at applied pressures > 7 GPa. The high- 
pressure phase obtained upon indentation is amorphous, as illustrated by the radial distribution 
functions calculated for different pairs of elements. In terms of elastic properties analysis, we have 
determined the elements of the stiffness tensor for a- and /?- eucryptite at the level of ReaxFF, and 
discussed the elastic anisotropy of these two polymorphs. Polycrystalline average properties of these 
eucryptite phases are also reported to serve as ReaxFF predictions of their elastic moduli (in the 
case of a-eucryptite), or as tests against values known from experiments or DFT calculations (/3- 
eucrypite). The ReaxFF potential reported here can also describe well single-species systems (e.g., 
Li-metal, Al-metal, and condensed phases of silicon), which makes it suitable for investigating struc- 
ture and properties of suboxides, atomic-scale mechanisms responsible for phase transformations, 
as well as oxidation-reduction reactions. 



I. INTRODUCTION 

Lithium Aluminum Silicate (LAS) glass ceramics have 
been investigated extensively over the last few decades 
owing to their exotic properties, such as small (or slightly 
negative) coefficient of thermal expansion and excep- 
tional thermal stability^— Such unique physical prop- 
erties make LAS ceramics suitable for a variety of ap- 
plications such as heat exchangers with high thermal 
shock resistance, high precision optical devices, telescope 
mirror blanks, and ring laser gyroscopes /3-eucryptite 
(LiAlSiO,i) is an important LAS glass ceramic material 
with a hexagonal crystal structure which can be viewed 
as a stuffed derivative of the high-temperature /3-quartz 
configuration^— The structure of /3-eucryptite confers 
it superionic conductivity of Li + ions along the c-axis, 
thus making it a potential electrolyte material for Li ion 
batteries^i^ 

From a more fundamental perspective, /3-eucryptite is 
known to undergo a reversible order-disorder transfor- 
mation at ~755 K which occurs via spatial disordering 
of the lithium atoms at high temperatures^ Recently, 
it was reported that /3-eucryptite undergoes a reversible 
pressure-induced transition to a metastable polymorph 
(called the e phase) at -0.8 GPaJiiS Apart from /3 
and the recently discovered high pressure e phases, there 
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are other known polymorphs of eucryptite, i.e. a 
and 7. Of these, a is the most stable under ambi- 
ent conditions and exists over a wide range of tempera- 
tures but is typically kinetically hindered.— On the other 
hand, 7-eucryptite is a metastable phase which coexists 
along with /3-eucryptite over a narrow range of temper- 
atures (1038-1103 K) and is, therefore, of lower practi- 
cal significance.— a-eucryptite has a rhombohedral crys- 
tal structure belonging to the R3 space group similar to 
phenakite and willemite; its thermodynamic, structural 
and physical properties have been studied . 21 ' 22 However, 
its elastic properties, in particular the single-crystal elas- 
tic constants are not yet known. Furthermore, the ex- 
tent of elastic anisotropy in a-eucryptite has also not yet 
been reported. Such a study could provide insight into 
the structure-property relationship in LAS glass ceram- 
ics and assist in designing composites with tailored elastic 
properties. 

The potential development work presented here has 
emerged from our long term goal of determining the 
atomic structure of the e phase and the atomic-scale 
mechanism responsible for the /3-to-e phase transforma- 
tion. While phase transformations can be directly evi- 
denced in DFT-based Carr-Parrinello molecular dynam- 
ics (CPMD) simulations (e.g., Refs. [H, Hi), the num- 
ber of atoms in the unit cell of /3 eucryptite makes it 
impractical to undertake such simulations in which usu- 
ally several unit cells are required along each spatial di- 
rection. We have therefore not resorted to CPMD ap- 
proaches, but turned to molecular dynamics (MD) sim- 
ulations based on empirical potentials. To describe the 
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interatomic forces acting during MD simulations, several 
empirical force fields (EFFs) have been proposed for LAS 
systems.— ~— These EFFs reproduce well short-range or- 
der, mechanical, and transport properties, but may not 
provide a sufficiently adequate description of phase trans- 
formations, medium range order, and vibrational den- 
sity of states. Recently, van Duin and co-workers 32,33 
developed a reactive force field (ReaxFF) based on a 
bond-order formalis m 34 i 35 in conjunction with a charge 
equilibration scheme—S. In this study, we have param- 
eterized ReaxFF for lithium aluminum silicates by fit- 
ting against formation energies, atomic configurations, 
and charge distributions of a number of representative 
clusters and equations of state of well-known condensed 
phases of oxides, silicates, and aluminates derived from 
DFT calculations. 

This article is organized as follows. Sec. UU de- 
scribes the methodology we adopted to parameterize the 
ReaxFF for lithium aluminum silicates, the DFT data 
used to construct the training set for the parametriza- 
tion of ReaxFF, and the details of these DFT calcula- 
tions. The parameters are given in a format compati- 
ble with the MD package LAMMPS^i Sec. HIl reports 
our results in for: heats of formation and geometric pa- 
rameters for a number of bulk phases; relative stability 
of eucryptite phases and amorphization of /3-eucryptite 
upon indentations; and elastic properties of the two most 
stable eucryptite polymorphs. We have determined the 
elements of the stiffness tensor for a- and /3-eucryptite at 
the level of ReaxFF, and discussed the elastic anisotropy 
of these two polymorphs. Polycrystalline average proper- 
ties of these eucryptite phases are also reported to serve 
as ReaxFF predictions of their elastic moduli (in the case 
of a-eucryptite) , or as tests against values known from 
experiments or DFT calculations (/3-eucrypite). Sec. |IV] 
summarizes the results and discusses the main successes 
and shortcomings of our ReaxFF parametrization for 
LAS systems. 



II. METHODOLOGY 

A. ReaxFF framework 

The energy contributions in ReaxFF are functions of 
bond orders, which allows for a better description of 
bond breaking and bond formation during simulations. 
Furthermore, Coulomb interactions are computed for 
every atom pair based on charges calculated at every 
time step using a charge equilibration scheme which al- 
lows it to describe covalent, metallic, and ionic systems 
equally well.— ^ The ReaxFF framework has been ap- 
plied successfully to predict the dynamics and reactive 
processes in hydrocarbons ) 32 ' 38 crack propagation in sil- 
icon crystals^ interfacial reactions in Si/Si-oxide^ 3 - and 
Al/Al-oxide,— surface reactions in ZnO, 41 oxygen-ion 
transport in Y-stabilized ZrC>2,— and phase transitions 
in ferroelectric BaTiC^. 43 The parameters of the ReaxFF 



potential are optimized by fitting against density func- 
tional theory (DFT) data for various bulk phases and 
atomic clusters. 

The formulation of ReaxFF is based on the concept 
of bond order ^ which describes the number of electrons 
shared between two atoms as a continuous function of 
their spacing. The bond order BO'^ associated with 
atoms i and j is calculated via 

BO' i4 = BOZ + BO'Z + BOT 
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where BO'^, BO'£ and BO'™ are the partial contribu- 
tions of ex, 7r- and double 7r-bonds involving the atoms i 
and j, Tij is the distance between i and j, Tq , rj, rj w are 
the bond radii of <r, n- and double 7r-bonds, respectively, 
and pbo are the bond order parameters. The bond orders 
BO'^ obtained from Eq. ([1]) are corrected to account for 
local overcoordination and residual 1-3 interactions by 
employing a scheme detailed in Ref. |45| . 

The total energy E of the system is expressed as 
the sum of partial energy contributions corresponding to 
bonded and unbonded interactions ! 32 ' 33 ' 38 
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where E bt ij is the energy of the i-j bond, E ov ^ and E un ^ 
are penalties for over- and under- coordination of atom 
i, E[ Pj i is the energy associated with lone-pair electrons 
around an atom i, E Vi ijk is the energy associated with 
the deviation of the angle subtended at j by atoms i and 
k from its equilibrium value, E v dw,ij and Ec,ij are the 
contributions from van der Waals and Coulomb interac- 
tions between i and j. 

The energy of the i-j bond is calculated using the cor- 
rected bond orders BOa as 



E b>ij = -D^BO^exp( Pbel (I - (BO^)p^-)) 
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where D° , and are the dissociation energies of a, 
it- and double 7r-bonds, while pbei,2 are the bond energy 
parameters. The contribution associated with lone pair 
electrons is calculated as: 
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where n/ Pi0p t is the optimum number of lone pairs for 
a given atom i and ni Pt i is the number of lone pairs 

around i calculated using the relation = ^-J + 
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where A„ e is the differ- 



ence between the number of outer shell electrons and the 
sum of bond orders around atom i and \_x\ is the greatest 
integer smaller than x. The penalty terms for overcoor- 
dination {E OVt i) and undercoordination (E un ^) of atom i 
can be written as 
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Vi is the valence of atom i, 



-j — T^lp.opt ^lp,J J 

Aj is the degree of overcoordination around the atom i 
which is corrected for the effect of broken electron pairs to 
obtain A l f c , and p's are over/under coordination param- 
eters. The energy contribution from the valence angles 
is written as: 
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is a shielding 



where fairy) = (i 
term included to avoid excessive repulsive interactions 
between bonded atoms and atoms containing a valence 
angle (1-3 interactions), is the depth of the potential 
well, r v dw is the van der Waal radius, p v dw and Jvdw 
are the van der Waals shielding parameters and T(ry ) is 
the Taper correction. The Coulomb interaction between 
atoms i and j is 



E c ,ij = T( rij )C- 



(8) 



where qi and qj are instantaneous charges of atoms i 
and j, C is the Coulomb constant and 7y is a shielding 
parameter included to avoid excessive repulsions due to 
overlap of orbitals at short distances. 

The atomic charges arc calculated at every iteration 
(time step) during minimization (MD) run using the 
Electronegativity Equalization Method (EEM)£ 6 This 
bond-order formalism coupled with the redistribution 
of charges through EEM enables the ReaxFF model to 
describe ionic, metallic, and covalent systems on equal 
footing,^H£8-ii.^i&-i2 W ith one unified (albeit com- 
plicated) formalism, ReaxFF has several advantages over 
other EFFs>^ 



E v , ijk = h iBO i:j ) f 7 [BO jk ) f 8 (A 3 -) F v (Q ijk ) (6a) 

F v i@ijk) =Pvi |l - cxp (-p V 2 (So - ©^jfe) 2 ) } (6b) 

where &ijk is the angle subtended at central atom j by 
the atoms i and fc, 0o is the equilibrium value for Qijk 
which depends on the sum of 7r-bond orders (i.e., BO™ 
and BO™™) around the atom j, fj and fs are functions of 
bond order and degree of overcoordination, respectively, 
and p v 's are valence angle parameters. 

All the terms on the right side of Eq. ([2]) except the van 
der Waals and Coulomb interactions depend on bond or- 
der through Eqs. ([J])-®. The bond orders are updated 
after every time step in a molecular dynamics simula- 
tion; such a formalism allows for a realistic simulation 
of dissociation and formation of bonds during a chemi- 
cal reaction and also provides a good description of the 
bulk phasesi 32 ' 33 ! 38 The pairwise non-bonded interaction 
terms, i.e., Coulomb and van der Waals interactions are 
evaluated for every atom pair irrespective of the geome- 
try and instantaneous connectivity. The van der Waals 
interaction of atoms i and j is evaluated as 



(i) The bond-order formalism provides a continuous 
description of formation and dissociation of bonds 
during a molecular dynamics simulation. 

(ii) Other interatomic potentials based on bond-order 
formalism like Tersoff 3 ^ and Brenner— do not ac- 
count for redistribution of charges. The EEM 36 
employed in ReaxFF allows the atomic charges to 
vary continuously with changes in coordination and 
bond order. 

(iii) The evaluation of individual contributions of a-, 
7T- and double n- bonds to the bond order allows 
ReaxFF to identify the hybridization state and the 
coordination of an atom based on the instantaneous 
geometry around that atom. 

(iv) The bond-order correction scheme enables ReaxFF 
to capture more accurately transition states dur- 
ing a reaction, provides a continuous transition 
between these intermediate states, and in turn 
describes reaction kinetics better than the other 
EFFs. 
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B. ReaxFF development 

The choice of specific partial energy contributions de- 
pends largely on the system of interest. For example, 
for ionic solids the angle bending and torsion terms have 
been set to zero] 43 ' 46 " — however, for covalent crystals 
these contributions cannot be neglected as shown, e.g., 
for the case of Si/Si02-— The partial energy contribu- 
tions used in the development of ReaxFF for Si/SiC>2 
system^ 3 - have been found adequate in the present study 
of Li/Al/Si/O as well [Eq. ©]. We optimized all ReaxFF 
parameters for the Li/Al/Si/O system by fitting against 
DFT-computed data. To ensure good transferability of 
the resulting Li/Al/Si/O parameters, we have included in 
the training set DFT-calculated data for a wide variety 
of well-known condensed phases and clusters, as listed 
below: 

(i) Equations of state {i.e. total energy versus volume) 
for pure Al (fee, hep, bec, sc and diamond) and 
for corundum (a-A^Oa), surface energy of the fee 
Al (111), charge distribution and dissociation en- 
ergies of a number of Al— O— H clusters; data from 
Ref. Mr 

(ii) Equations of state of Li (bee, fee, hep, diamond, 
sc), LiH with sodium-chloride structure, dissocia- 
tion energies and charge distributions in Li2, LiH 
and LiH2 clusters; data from Ref. I4IH 

(iii) Equations of state of Si (sc, diamond, /3-Sn), 
SiC>2 (a-quartz, trydimite, coesite, a-crystobalite, 
stishovite), dissociation energies of single and dou- 
ble bonds of Si— Si and Si— O in Si/O/H clusters, 
energies of various Si/O/H clusters as a function of 
valence angles Si— O— Si, O— Si— O and Si— Si— Si 
and distortion energies of rings of Si/O/H clusters; 
data from Ref. [H. 

(iv) Equations of state of Li-silicates: (a) Li 2 Si03 (or- 
thorhombic) (b) Stable Li 2 Si20s (monoclinic) and 
(c) Metastable Li2Si2 05 (orthorhombic); data from 
Ref. [M 

In addition to using DFT data sets from earlier works, 
we calculated the equations of state of the following con- 
densed phases within the framework of DFT using the 
computational details listed in Sec. Ill CI 

(v) Li-oxides: a-Li20 (cubic) 51 and Li202 
(hexagonal) £2- 

(vi) Li-aluminates: Three polymorphs of LiA102 
namely, (a) a (rhombohedral) 53 , (b) /3 
(orthorhombic) ^ and (c) 7 (tetragonal) j 5 - 5 - 

(vii) Al-silicates: Three polymorphs of A^SiOs 
namely, (a) Andalusite (orthorhombic) t 56 i 57 (b) 
Sillimanite (orthorhombic), 5 ?' 58 and (c) Kyanite 
(triclinic)^^ 



In order to account for anisotropy, computational su- 
percells of these phases were subjected to different types 
of strain (depending on the crystal symmetry) when com- 
puting their energy as a function of cell volume. The 
lattice of a crystal is described by three lattice vectors an 
(i = 1, 2, 3) whose magnitudes are the lattice parameters 
di. The cubic phases (a\ = 02 = 03) were strained triaxi- 
ally, i.e., all the three lattice vectors (ai) were all equally 
strained. The tetragonal and hexagonal phases (ai = 
^ 03) were deformed by two types of strains, namely (a) 
biaxial: a\ and 0,2 were strained simultaneously by the 
same amount while keeping 03 fixed at its experimen- 
tal value, and (b) uniaxial: 03 was strained while keep- 
ing a\ and 02 fixed. The lattice vectors of phases with 
orthorhombic symmetry or lower [a\ /«2/ 03) were 
strained individually keeping the other two unstrained, 
which leads to three distinct uniaxial strains correspond- 
ing to three lattice vector directions a$. In all the cases, 
the limits of strains range from —40% (compressive) to 
+20% (tensile). 

For all the phases listed above [items (i)-(vii)], we com- 
puted the heats of formation AHf as functions of volume 
for the different types of strains. The heat of formation of 
a general compound of unit formula (u.f.) LifcAliSi m O„ 
(fc, Z, m, n integers > 0) at a volume V for a particular 
type and value of strain can be evaluated from DFT to- 
tal energy calculations as: 

AH f (V,e) = E Uk Ai lS i m o n (V,e)-kE u 

ft 

-IE ai - mE Si - -E . 2 (9) 

where -ELi fc Ai(Si m O„ is the total energy of a given volume 
V of the phase LifcAl;Si m O„ subjected to a particular 
strain e. The energies of the constituent elements Li, Al, 
Si and O in Eq. ((9]), i.e., Eu, Eai, E$i, and Eo 2 , are 
those of the most stable phases at equilibrium calculated 
by DFT. 

The training set data were used to parameterize the 
ReaxFF using the successive one-parameter search tech- 
nique described by van Duin et al£Q These parameters 
are tabulated in Appendix^ and are also made available 
as a data file.— 



C. Details of the DFT calculations 

The computational supercell for each phase in the 
training set described in Sec. Ill Bl consisted of one prim- 
itive unit cell. The total energy DFT calculations 
were performed within the framework of the general- 
ized gradient approximation (GGA), using the projector- 
augmented wave (PAW) formalism 6 -^ as implemented in 
the ab-initio simulation package VASPi 63 ' 64 The atomic 
coordinates were relaxed using a conjugate gradient al- 
gorithm until the force components on any atom were 
smaller than 0.01 eV/A. The exchange-correlation was 
described by the Perdew-Wang functional, 65 which is a 
typical choice for ceramics oxide systems (e.g., Ref. IHol 
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FIG. 1: (Color online) Equations of state of various phases of (a, b) Li oxides, (c, d) Li aluminates, (e, f) Al silicates, and 
h) Li silicates as calculated using DFT [panels (a), (c), (e), (g)] and ReaxFF [panels (b), (d), (f), (h)]. 



I55T ). The plane wave energy cutoff was set to 500 
eV, which performs satisfactorily for similar ceramic 
systems^ The Brillouin (BZ) zone was sampled with 
a T-centered Monkhorst-Pack grid. For the oxides and 
aluminates of lithium, we used 8x8x8 fc-point grids 
which amount to 1024 irreducible fc-points for oxides and 
256 fc-points for aluminates. A 4 x 4 x 4 fc-point grid 
was found sufficient for the aluminum silicate phases (32 
irreducible fc-points), and a 3 x 3 x 3 fc-point grid was se- 
lected for the eucryptite phases (14 irreducible fc-points). 



These grids were chosen on the basis of convergence tests 
conducted for different BZ samplings for different phases. 

III. RESULTS 

A. Heats of formation 

The set of parameters obtained by the technique de- 
scribed in Sec. II were validated by comparing the struc- 
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TABLE I: Heats of formation at equilibrium (AHf) of se- 
lected phases calculated using DFT and ReaxFF at K. For 
comparison, experimental values at 298 K are also provided 
wherever available. 



Phase 




AH° f (kcal mol" 1 ) 




DFT 


ReaxFF 


Exp 


a-Li 2 


-147.67 


-145.82 


-143.10° 


7-L1AIO2 


-292.19 


-299.29 


-284.37 b 


Andalusitc 


-635.90 


-635.35 


-619.42 c 


Li2Si03 


-406.04 


-387.31 


-395.77 d 


a-LiAlSi0 4 


-529.25 


-526.48 


-512.53 e 


/3-LiAlSi0 4 


-528.01 


-525.51 


-506.18 / 


7-LiAlSi0 4 


-524.93 


-517.66 





a Ref.[6(J ^Ref.^; c Ref.[6g; d Ref.[69|; e Ref.|2j; -fRef-IZOl 



tures and the heats of formations for various phases cal- 
culated by ReaxFF with those known from experiments 
or from DFT calculations. As a preliminary test, the 
heats of formation as functions of volume for the vari- 
ous phases used in the training set calculated by ReaxFF 
were compared in Fig. [T] with their DFT counterparts. 
Figure Q] shows a generally good qualitative agreement 
between the ReaxFF and the DFT curves in terms of 
equilibrium volumes and the relative phase stabilities at 
these volumes. Furthermore, Table U shows that the 
ReaxFF heat of formation results are also in good agree- 
ment with experimental data on selected oxides, alumi- 
nates, and silicates at equilibrium. However, in the defor- 
mation regimes lying outside equilibrium (particularly in 
tension) the energetic ordering of Li oxides [Fig. [IJa,b)], 
Li aluminates [Fig. DJc,d)], and Al silicates [Fig. [TJe,f)] 
at the ReaxFF level does not preserve so well the DFT 
ordering. This is most likely due to the choice of deforma- 
tion range (Sec. Ill Bj) . which contains more data points 
in compression than in tension. The following subsec- 
tions contain more tests of the performance of ReaxFF 
concerning the structure, stability, and elastic properties 
of LAS ceramics. 



B. Structural parameters 

Table iHl compares the lattice parameters for a number 
of selected phases calculated using ReaxFF with those 
from DFT calculations and from experiments. These 
lattice constants were calculated by optimizing the com- 
putational supercell of each phase with respect to all 
independent lattice parameters that describe its crystal 
structure. Table HI1 shows that the lattice parameters cal- 
culated using ReaxFF are all within ^5% of the values 
reported in literature using DFT or experiments. In or- 
der to establish that the structures of bulk phases are 
faithfully described by ReaxFF, we have also checked 
the independent fractional coordinates of the atoms in 
the optimized supercells. For example, Tables Hill and HVl 



show these fractional coordinates for a- and /3-eucryptite, 
respectively. As shown in these tables, the agreement 
between the ReaxFF-predicted values for fractional co- 
ordinates and the experimental ones is very good, which 
illustrates that ReaxFF predicts the structure of bulk 
phases accurately. 



C. Stability of eucryptite phases 

There are three well-known crystalline polymorphs of 
eucryptite, a, (3, 7; of these, a is the most stable phase 
under ambient conditions but is kinetically hindered. 19 
Fig. [5] shows the equations of state (Energy vs Volume 
curves) for these polymorphs calculated using ReaxFF 
and DFT at K. The minimum of a calculated energy vs 
volume curve for a given phase represents its equilibrium 
state. For convenience, all the energies reported in Fig. [2] 
are given relative to the energy of the most stable phase 
at its equilibrium volume for DFT and ReaxFF. Fig. [2] 
shows that ReaxFF predicts the same order of stability 
for the three polymorphs of eucryptite as do our DFT 
calculations. This order is consistent with experimental 
observations i 3 ' 11 ! 17 ' 21 
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FIG. 2: (Color online) Equations of state of various phases of 
eucryptite calculated using (a) DFT and (b) ReaxFF. Both 
the techniques predict the same order of stability of the three 
polymorphs with a being the most stable phase in each case. 



/3-eucryptite, the most technologically relevant of the 
three polymorphs, has a open structure which collapses 
at sufficiently high applied pressures . 1 ' 2 ' 18 Recently, it 
was observed that /3-eucryptite begins to amorphize at 
pressures above ~5 GPa^S To test the ability of ReaxFF 
to capture phase transitions, we have studied the evolu- 



7 



TABLE II: Comparison of calculated lattice parameters of selected phases using ReaxFF with those available in literature 
determined from DFT calculations and from experiments. The numbers of unit formulae per unit cell are provided in brackets. 



Phase Structure Space Group Formula Lattice DFT ReaxFF Exp 

parameter (A) 





Cubic 


nun ^TTJ 


T,iiO (41 




4.631° 


4.738 


4.622 6 


V— I 11 n l\ 1 2 




P4 1 2 1 2 






5.223 c 


o.oou 


^ 1 f>Q d 










c 


6.309 


6.234 


6.ZD8 


Andalusite 


Orthorhombic 


Pnnm 


AbSiOs (4) 


a 


7.753 e 


7.632 


7.798 ; 








b 


7.844 e 


7.916 


7.903 / 










c 


5.477 e 


5.727 


5.557 J 


Li2Si03 


Orthorhombic 


Cmc2 1 


Li 2 Si0 3 (4) 


a 


9.487 9 


9.335 


9.392 ft 










b 


5.450 9 


5.431 


5.397 h 










c 


4.713 9 


4.861 


4.660' 1 


a-eucryptite 


Trigonal 


R3 


LiAlSi0 4 (18) 


a 


13.656 


13.448 


13.532' 










c 


9.158 


8.981 


9.044 i 


/3-eucryptite 


Hexagonal 


P6 4 22 


LiAlSiQ 4 (12) 


a 


10.594 J 


10.568 


10.497 fc 










c 


11.388 J 


11.763 


11.200 fc 



a Rcf.[5l|; b Ref.l7l|; c Ref.l55|; d Ref.l7^; e Ref.l7^; / Ref.|5^; 9 Ref.l7^; ''Ref-lZfil; ^Ref.^; JRef.^; fe Ref.lll| 



TABLE III: Fractional coordinates of atoms in a unit cell of a- 
eucryptite calculated using ReaxFF at K. The experimental 
values from Ref. [2l] at 298 K are provided for comparison. 



TABLE IV: Fractional coordinates of atoms in a unit cell of /?- 
eucryptite calculated using ReaxFF at K. The experimental 
values from Ref. 3 at 298 K arc provided for comparison. 







ReaxFF 




Experiment 


Atom 


X 


y 


z 


X 


y 


z 


Li(l) 


-0.016 


-0.806 


-0.752 


-0.017 


-0.811 


-0.749 


Li(2) 


0.022 


0.814 


0.749 


0.021 


0.812 


0.754 


Si(l) 


0.531 


0.876 


0.753 


0.530 


0.880 


0.750 


Si(2) 


0.876 


0.348 


0.918 


0.876 


0.344 


0.916 


Al(l) 


-0.533 


-0.883 


0.754 


-0.530 


-0.882 


-0.749 


Al(2) 


-0.878 


-0.342 


-0.914 


-0.875 


-0.345 


-0.916 


0(1) 


-0.748 


-0.208 


-0.897 


-0.753 


-0.210 


-0.890 


0(2) 


0.764 


0.211 


0.903 


0.766 


0.216 


0.898 


0(3) 


-0.741 


-0.202 


-0.594 


-0.733 


-0.199 


-0.593 


0(4) 


0.734 


0.198 


0.571 


0.733 


0.199 


0.576 


0(5) 


-0.097 


-0.886 


-0.931 


-0.105 


-0.888 


-0.937 


0(6) 


0.090 


0.879 


0.947 


0.096 


0.881 


0.946 


0(7) 


-0.669 


-0.009 


-0.751 


-0.664 


-0.009 


-0.749 


0(8) 


0.656 


-0.004 


0.753 


0.655 


-0.004 


0.750 



tion of /3-eucryptite under a rigid spherical indenter us- 
ing MD simulations^ An orthorhombic simulation box 
of dimensions 41.99A x 72.73A x 56A containing 13440 
atoms was used to simulate the crystal, which was in- 
dented down the z axis (i.e., the [OUT] crystal direction). 
Periodic boundary conditions were applied in the direc- 
tions perpendicular to the indentation force. The atoms 
that have z coordinates within 12A of the lowest z value 
(of all atoms) were kept fixed during MD runs in order 
to simulate the underlying bulk. The initial structure 
was relaxed at K, and then thermalized at 300 K for 
30 ps; the time-step used in the MD runs was 1 fs. Af- 
ter thermalization, the top face was indented at a rate of 
0.065A/psby a rigid spherical indenter of radius R = 14 A 



Atom 






ReaxFF 








Experiment 






X 


y 




z 




X 


y 




z 


Li(l) 





.000 


0.000 





.500 


0. 


.000 


0.000 





.500 


Li(2) 





.500 


0.000 


0. 


.000 


0. 


.500 


0.000 





.000 


Li(3) 





.500 


0.000 





.327 


0. 


.500 


0.000 





.328 


Si(l) 





.248 


0.000 





.000 


0. 


.248 


0.000 





.000 


Si(2) 





.251 


0.502 


0. 


.000 


0. 


.247 


0.494 





.000 


Al(l) 





.258 


0.000 





.500 


0. 


.250 


0.000 





.500 


Al(2) 





.250 


0.499 





.500 


0. 


251 


0.501 





.500 


0(1) 





.115 


0.201 


0. 


.248 


0. 


112 


0.199 





.242 


0(2) 





.101 


0.696 





.263 


0. 


097 


0.699 





.259 


0(3) 





.604 


0.704 


0. 


.262 


0. 


.597 


0.705 





.264 


0(4) 





.605 


0.202 





.258 


0. 


.608 


0.201 





.249 



which applies a radial force Fi on atom i given by: 

f, = < ,:. (10) 



~k{n - R) 2 iin<R 
if n > R 



where k is a force constant (k — 76.32 kcal/mol A 3 ), and 
Ti is the distance between the center of the atom i and 
that of the indenter. 

During indentation simulations, we have not found any 
new phase at indent pressures smaller than 7 GPa even 
though the e phase has been reported^ to occur at ~0.8 
GPa. One reason for which we do not observe the e phase 
in these simulations is that the pressure is not applied 
hydrostatically (as it was in experiments 17 ), and such 
anisotropic application of external pressure may lead to 
different phase transitions, 2 - i.e., different onset pressure 
or different phases. We have observed the formation of 
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a denser but disordered phase in the vicinity of the in- 
denter. Zhang et al. 2 ® carried out in-situ X-ray experi- 
ments on polycrystalline samples of /3 eucryptite and re- 
ported that amorphization begins at a pressure ^5 GPa 
and completes at a pressure of 17 GPa. Our indenta- 
tion MD simulations predict a higher onset pressure for 
amorphization, ~7 GPa. However, it should be noted 
that the simulations were carried out on single crystal /3 
eucryptite; this eliminates the defects, porosity, or grain 
boundaries from our starting phase which would have 
acted as nucleation sites for the formation of the amor- 
phous phase. Consequently, in the case of our simula- 
tions, there is an increased barrier towards amorphiza- 
tion, which is reflected in the increased onset pressure. 
At the ReaxFF onset pressure of 7 GPa, only regions 
near the indent amorphize, while those far away from it 
remain crystalline. As the indent pressure is increased, 
the amorphizcd region grows and there is a range of pres- 
sures over which amorphization proceeds: this finding is 
similar to what of Zhang et al. have found for polycrys- 
talline samples^ 

We have analyzed the amorphous phase by studying 
radial distribution functions (RDF) for pairs of differ- 
ent types of atoms in the disordered region. Specifically, 
the RDFs <7a— b were evaluated for Si-O, Al-O, Li-0 
and Li-Li pairs after indentation proceeded to different 
depths h. Fig. [3] shows RDFs evaluated prior to the in- 
dentation, compared to those calculated after indenta- 
tion to h = 12 A; at this depth, we evaluated the con- 
tact pressure at ~10 GPa. The RDFs for all the type 
pairs considered show, prior to indentation, well-defined 
peaks at characteristic distances of /3-eucryptite. Un- 
der pressure, the first peak (r = 1.6 A) of gsi-o( r ) [see 
Fig. [3ta)] broadens somewhat and decreases in intensity 
compared to that of the crystalline j3 phase. The peak 
corresponding to the second-nearest neighbor (r = 4.1 A) 
is very broad, while peaks at higher distances are not de- 
fined. A similar behavior of the RDF was observed in the 
amorphous phase obtained from high pressure MD sim- 
ulations of /3-crystobalite (Si02).~ Fig. (3^b) shows the 
RDF for Al-O, which also exhibits the tell-tale signs of 
a disordered phase under pressure; the broadening of the 
first peak, along with the disappearance of the higher- 
order peaks, has also been observed during amorphiza- 
tion of SiG^ and a-quartz^ Significant changes in the 
RDFs occur at contact pressures >7 GPa and indicate 
amorphization, which is also apparent from direct visu- 
alization of the structure. The pressure necessary for 
the onset of amorphization is consistent with empirical 
observations^ 

Interestingly, an additional feature is exhibited by the 
RDFs of Li-0 and Li-Li pairs [Fig. G2c,d)]. The first 
peaks for the Li-Li and Li-0 pairs are shifted signifi- 
cantly to lower distances in the high pressure phase as 
compared to the initial crystal [Fig. G3c,d)]. For the Li- 
O pairs, the first peak shifts from 2 A to 1.67A under 
pressure [Fig. Etc)], which is close to the typical Li-0 
bond length of 1.606A. 80 Fig. [2Id) shows that the small- 



TABLE V: Predicted stiffness constants CV, (in GPa) of a- 
eucryptite using ReaxFF at K. 



Cii 


Cl2 


Cl3 


C24 


Cu 


C33 


C44 


ReaxFF 131.86 


67.92 


25.18 


1.96 


1.28 


175.24 


37.29 



est most probable Li-Li spacing at high pressure is ~ 
2.88A which is closer to the experimental value of the 
bond length (3.04A) in Li-metal^ than it is to the low- 
est Li-Li distance (3.8A) in /3-eucryptite. This suggests 
the existence of Li-Li bonds in the high pressure phase, 
which were not present in the crystalline phase; by check- 
ing the atomic structure details of the amorphized phase, 
we have confirmed the presence of direct Li-Li bonds 
and have also found bonds in which two Li atoms are 
"bridged" by on O atom in a triangular configuration. 
These newly formed, compressed Li-Li bonds and the 
shortened Li-0 bonds formed under pressure suggest den- 
sification in the amorphous phase. 



D. Elastic properties of eucryptite phases 

To assess the performance of ReaxFF in predicting 
elastic properties, we computed the elements of elastic 
stiffness tensor Cy for two polymorphs of eucryptite, 
a and j3, by employing the technique outlined in Ap- 
pendix [B] We have found that the stiffness tensors for 
both eucryptite phases are positive definite, which means 
that at the ReaxFF level the Born stability criterion^ 
is met. The seven independent elastic constants of a- 
eucryptite (rhombohedral structure) were calculated us- 
ing ReaxFF at K and are listed in Table fVl to the best 
of our knowledge, so far there are no reports of elastic 
constants in the literature for this phase. 

We have also calculated the five independent elas- 
tic constants of /3-eucryptite (hexagonal structure) at 
K predicted by ReaxFF and listed them in Table I VII 
Haussiihl et al£2- have measured these elastic constants 
using an ultrasonic technique at ambient temperature, 
293 K. To compare the elastic constants predicted by 
ReaxFF in our study and those predicted by DFT (also 
at K)2l with the experimental values, we have extrapo- 
lated the measured values of CV,- to K using thermoelas- 
tic constants = (flogCy / dT M As shown in Table IVTl 
the ReaxFF elastic constants are in good agreement with 
experiment; with the exception of C12, all the calculated 
constants are within ~ 30% of experimental values ex- 
trapolated to K. These ReaxFF values are also con- 
sistent with those predicted by DFT and reported in an 
earlier workiSi 

The a and j3 polymorphs of eucryptite are known 
to possess highly anisotropic physical properties. The 
overall elastic anisotropy of hexagonal and rhombohedral 
crystals is usually assessed through three ratios, C11/C33, 
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4 
3 
2 
1 


4 
3 
2 
1 


0.6 
L 0.4 
0.2 





(a) 



amorphous - 

Si-0 






6 

r(A) 



10 



12 



FIG. 3: (Color online) Pair distribution functions (<?a-bM) 
for (a) Si-O, (b) Al-O, (c) Li-O and (d) Li-Li pairs in /3- 
eucryptite (black lines) and in the phase obtained under a 
spherical indent (red lines) at an applied contact pressure 
~ 10 GPa. The broadening of the peaks corresponding to 
higher order neighbors and lowering of the nearest-neighbor 
distances in Li-O and Li-Li pairs at high pressures indicates 
that the new phase formed under the indent is amorphous. 



C12/C13 and 2644/(611 — C12), whose deviations from 
unity serve as measures of the anisotropy in the crystals 
being studied. Table I VIII lists the anisotropy ratios for 
rhombohedral a and hexagonal /3-eucryptite using the 
Cij predicted by ReaxFF. 

The anisotropy of eucryptite polymorphs manifests, 
expectedly, in the Young's modulus E as well as in 
other elastic properties. The direction-dependence of 
Young's modulus can be derived from the elastic con- 
stants, and we show it here as a way to directly visualize 
the anisotropic character of the Young's modulus (Fig.0|. 
The Young's modulus for a rhombohedral crystal in the 



TABLE VI: Comparison of the calculated stiffness constants 
dj (in GPa) of 3-eucryptite at K with the experimental 
data from Ref.|83 extrapolated to K using the thermoelastic 
constants Ty = d log CijjdT. The uncertainty in any of the 
experimental values (Exp) is smaller than 2.5 GPa. 





C11 


C12 


C13 


C33 


C44 


DFT a 


165.64 


70.98 


78.59 


132.83 


58.68 


ReaxFF 


178.92 


102.77 


118.28 


181.26 


47.37 


Exp 


176.3 


68.5 


89.8 


139.9 


61.2 


Tij (W- 3 /K) 


-0.14 


0.13 


-0.27 


-0.42 


-0.24 



a Ref.lM K 



TABLE VII: Anisotropic factor ratios for a and j3 eucryp- 
tite (LiAlSi04) evaluated using single crystal elastic constants 
(Cij) predicted by ReaxFF in the present study. For /?- 
eucryptite, these ratios are also calculated using dj known 
by DFT and experiments for comparison. 



Phase 


Technique 


C11 

C33 


C12 
C13 


2C44 
Cn — C12 


a-LiAlSi0 4 


ReaxFF 


0.7525 


2.6974 


1.1664 


/3-LiAlSi0 4 


DFT a 
ReaxFF 
Exp 6 


1.2470 
0.9871 
1.2602 


0.9032 
0.8688 
0.7628 


1.2398 
1.2441 
1.1354 



6 From dj in Ref. 



extrapolated to K. 



R3 space group along a crystallographic direction of di- 
rection cosines l\, I2, h can be expressed in terms of 
elastic compliance constants as£^ 

i = (l-2 3 2 ) 2 S 11 +Z 3 4 S3 3 + ^(l-^)(2Si3 + S44) 

+2Z 2 Z 3 (3Z 2 - %)S U + 2hl 3 (3ll - ll)S 25 , (11) 

where Sij are the elements of elastic compliance matrix, 
S given by the inverse of the elastic stiffness matrix, C 
i.e., S = C _1 . For hexagonal crystals, the directional 
dependence of E is given by^ 



1 

E 



(1 - ijfSn + liS 33 + Z 2 (l - l 2 3 )(2S 13 + S 44 ) (12) 



Figure SJa) shows the variation of Young's modulus 
of ev-eucryptite with the angle 9 between a given crys- 
tallographic direction and the z-axis, for three different 
planes containing the z-axis; these planes are yz (h = 0), 
xz (I2 = 0) and the plane containing the first bisector of 
the xy plane, h = h = sin 8/ \p2. The three polar plots in 
Fig. |4ja) were generated using Eq. (fTTj) and the ReaxFF 
elastic constants in Table [V] 

The Young's modulus of /3-eucryptite depends only on 
the angle between a given direction and the z-axis 
(crystallographic c-axis), owing to the symmetry of a 
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FIG. 4: (Color online) Polar plots illustrating the directional 
dependence of Young's Modulus E for (a) a-eucryptite us- 
ing the elastic constants predicted by ReaxFF in three dif- 
ferent crystallographic planes containing the z-axis and (b) 
/3-eucryptite using the elastic constants predicted by ReaxFF 
and those known by DFT and experiments. 



hexagonal crystal. Figure [4] shows the dependence of 
Young's modulus of /3-eucryptite on calculated using 
the ReaxFF elastic constants and those known from DFT 
and experiments (Table [VT]). The variation with 8 of the 
Young's modulus of /3-eucryptite calculated using elastic 
constants predicted by ReaxFF follows the same trends 
as the E calculated using the elastic constants known by 
DFT or experiments. We now focus on elastic properties 
corresponding to polycrystalline eucryptite phases. 

The theoretical average bulk (B) and shear (G) elas- 
tic moduli of polycrystalline a and /3-eucryptite can 
be derived from their single-crystal elastic constants. 
There are two well-known approximations typically used 
to evaluate the polycrystalline elastic moduli, namely 
the Voigti^ and Reuss 86 methods, which provide upper 
bounds (identified by the subscript V) and lower bounds 
(subscript R), respectively, for the bulk and shear mod- 
uli. For rhombohedral and hexagonal crystal systems, 
the polycrystalline bulk and shear moduli can be ex- 
pressed in terms of the single-crystal elastic constants 



TABLE VIII: Average values of bulk moduli (B, in GPa) and 
shear moduli (G, in GPa) using the Voigt, Ruess, and Hill's 
approximations for polycrystalline eucryptite phases derived 
from their single- crystal elastic constants dj. The Young's 
moduli (Epdy, in GPa) and Poisson ratios (v po i y ) are evalu- 
ated using Eqs. (fl7)) . 





Q-LiAlSi0 4 




^-LiAlSi0 4 




ReaxFF 


DFT" 


ReaxFF 


Exp 


B v 


75.06 


102.27 


135.31 


109.85 


B R 


75.06 


101.55 


134.88 


109.55 


B h 


75.06 


101.91 


135.09 


109.70 


Gy 


42.69 


48.67 


39.88 


51.55 


Gr 


38.44 


46.08 


38.48 


47.10 


G h 


40.56 


47.37 


39.18 


49.32 


Epoly 


103.12 


123.05 


107.18 


128.69 




0.27 


0.30 


0.37 


0.31 



"From dj in Ref. 
°From dj in Ref. | 



as 



extrapolated to K. 



— = 2S n + 2S 12 + 45i 3 + S-. 







By 


"5< 


1 




= 25 


B^. 






1 


Gy 




30 


1 


1 


Gr 


15 



:i.'S 



(13) 
(14) 



The Hill values (Bh, Gh) of the bulk and shear moduli 
are the arithmetic averages of the corresponding Voigt 
and Ruess bounds, and are considered the best estimates 
of these polycrystalline moduli.— The polycrystalline 
Young's modulus E po i y and Poisson's ratio (v po i v ) can 
be obtained through the relations applicable to isotropic 
materials,— 



Epoly — 



9B H G H 
3Bh + Gh 



3Bh — 2Gg 
2{3B H + G H )' 



(17) 



Table I Villi lists the average elastic moduli of polycrys- 
talline a and (3 eucryptite derived from the single-crystal 
constants Gy (Tables [V] and IVI|) through the relation- 
ships in Eqs. (|T3l) - (fT7l) . For a-eucryptite, the polycrys- 
talline bulk modulus evaluated using the single crystal 
elastic constants predicted by ReaxFF (75.06 GPa) is 
in excellent agreement with an earlier measurement (74 
GPa) of Fasshauer et al?" 1 Furthermore, the polycrys- 
talline elastic constants of /3-eucryptite calculated from 
the single crystal elastic constants predicted by ReaxFF 
is in good agreement with those calculated using the elas- 
tic constants known by DFT and experiments (refer to 
Table ELLE}. 



11 



IV. SUMMARY AND CONCLUSION 

In conclusion, we have developed a reactive force field 
for lithium aluminum silicates and used it to describe (i) 
the atomic structure and heats of formation of several 
oxides, silicates and aluminates, (ii) the relative stabil- 
ity of three crystalline eucryptite polymorphs and the 
response of /3-eucryptite under indentation, and (iii) the 
anisotropic and polycrystalline-averaged elastic proper- 
ties of eucryptite phases. 

Successes. We have found that structural proper- 
ties and heats of formation for selected condensed phases 
agree well with the results of DFT calculations and with 
experimental reports. In terms of applications to the sta- 
bility of eucryptite phases, we have verified that the order 
of the stability of three well-known polymorphs predicted 
by ReaxFF is the same as that obtained from DFT calcu- 
lations and that known from experiments. The response 
of /3-eucryptitc to pressure is the formation of a denser 
and disordered phase which we characterized by a set of 
radial distribution functions and comparisons with con- 
densed phases. In terms of elastic properties analysis, we 
have determined the elements of the stiffness tensor for a- 
and /3- eucryptite at the level of ReaxFF, and discussed 
the elastic anisotropy of these two polymorphs. Polycrys- 
talline average properties of these eucryptite phases are 
also reported to serve as ReaxFF predictions of their elas- 
tic moduli (in the case of a-eucryptite), or as tests against 
values known from experiments or DFT calculations (/3- 
eucrypite). In addition to the elaborate but physically- 
motivated description of the bond order formalism cou- 
pled with the EEM scheme, the novel aspects/results of 
this work include the ability of ReaxFF to predict the for- 
mation of an amorphous phase under pressures exceeding 
7 GPa, and the prediction of all elastic properties of a- 
eucryptite -which is the most stable LiAlSi04 phase at 
room temperature and ambient pressure. 

Shortcomings. We noted in Sec. IIIA that in the de- 
formation regimes far outside equilibrium, the ReaxFF- 
predicted order of phase stability may not match the 
DFT predictions, especially in the tensile regimes. This 
problem is likely to manifest during reaction calculations 
at the level of several atoms, molecules, or small clusters, 
but may not easily manifest in large-scale MD simula- 
tions because fracture in tensile regimes will probably 
occur before any phase transformation. The values of 
the ReaxFF elastic constants Cy for /3-eucryptite com- 
pare reasonably well with those predicted by other em- 
pirical force fieldsi^Zr— These values are not of superior 
accuracy, as they deviate by about 30% from the ex- 
perimental values. Still, the values of predicted by 
ReaxFF ( Table IVTj) deviate from experiments by amounts 
that are very similar to the deviations made by the Pe- 
done force field in predicting the elastic constants for spo- 
dumenc (LiAlSi20e)>^ However, we found that PFF— 
and a core-shell model potential developed by Winkler et 
al. (THB) 28 describe the elastic properties of aluminum 
silicate phases (especially those of andalusite, A^SiO,^) 



much better than ReaxFF. While this observation seems 
to place ReaxFF at a disadvantage, it should be noted 
that the PFF and other models were obtained by fit- 
ting against experimental values of the elastic properties 
of binary oxides and silicates; 28 ' 31 while these properties 
were not a part of the training set used to determine the 
ReaxFF parameters in our study. 

Concluding Remark. The ReaxFF potential re- 
ported here can also describe well single-species systems 
(e.g., Li-metal, Al-metal, and condensed phases of sili- 
con), which makes it suitable for investigating structure 
and properties of suboxides, atomic-scale mechanisms re- 
sponsible for phase transformations, as well as oxidation- 
reduction reactions. Based on the results of indentation 
on /3-eucryptite and elastic properties of a-eucryptite 
reported here, we believe that the parametrization of 
ReaxFF for Li-Al silicates will help provide fundamental 
understanding of other interesting phenomena in LAS 
glass ceramics, especially in regard to the atomic scale 
mechanisms underlying the pressure induced /3-to-e phase 
transformation where direct dynamic simulations at the 
level of DFT are currently intractable. 
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Appendix A: ReaxFF parameters for Li-Al-Si-O 

systems 

The ReaxFF parameters for the Li-Al-Si-0 systems 
determined in the present study are listed in Ta- 
bles [O-EvIIIl 

TABLE A.I: General Parameters 



TABLE A. IV: Van der Waals interaction parameters. 



Atom r vdw (A) Djj (kcal/mol) a -y vdw (A) 



Li 


1.6121 


0.2459 


10.8333 


1.4649 


Al 


2.3738 


0.2328 


9.4002 


1.6831 


Si 


1.8951 


0.1737 


11.3429 


5.2054 





2.3890 


0.1000 


9.7300 


13.8449 



Parameter 


Value 


Description 


Pbocl 


50.0000 


Bond order correction 


Vbocl 


9.5469 


Bond order correction 


P3 


50.0000 


Overcoordination 


Pi 


0.6991 


Overcoordination 


P6 


1.0588 


Undercoordination 


P7 


12.1176 


Undercoordination 


P8 


13.3056 


Undercoordination 


Plpl 


6.0891 


Lone pair parameter 


Pv7 


33.8667 


Valence undercoordination 


PvS 


1.8512 


Valence angle 


Pv9 


1.0563 


Valence angle 


PvW 


2.0384 


Valence angle 


PvdWl 


1.5591 


van der Waals shielding 


BOcut 


0.0010 


Bond order cut-off 



TABLE A.V: Bond parameters. The bond dissociation ener- 
gies D1 and DJ* are in kcal/mol while piei, Pbe2 and pi 
are unitless 









n 1 


JJTV7V 


Pbel 


Pbe2 


Pi 


o- 


O 


142.2858 


145.0000 


50.8293 


0.2506 


-0.1055 


0.3451 


Si- 


-0 


274.8339 


5.0000 


0.0000 


-0.5884 


-0.2572 


9.9772 


Si- 


-Si 


70.9120 


54.0531 


30.0000 


0.4931 


-0.8055 


0.2476 


Al- 


-0 


181.1998 


0.0000 


0.0000 


-0.2276 


-0.3500 


0.2086 


Al- 


-Si 


0.0000 


0.0000 


0.0000 


1.0000 


0.0000 


0.5000 


Al- 


-Al 


34.0777 


0.0000 


0.0000 


0.4832 


-0.4197 


6.4631 


Li- 


-o 


78.3666 


-0.0200 


0.0000 


-1.0000 


-0.2500 


0.2022 


Li- 


-Si 


0.0000 


0.0000 


0.0000 


1.0000 


0.0000 


0.5000 


Al- 


-Li 


0.0000 


0.0000 


0.0000 


1.0000 


0.0000 


0.5000 


Li- 


-Li 


42.9780 


0.0000 


0.0000 


0.3228 


0.0000 


1.7161 



TABLE A. II: Atom parameters. All the parameters except 
Pip2 (kcal/mol) are unitless 



Atom 




V| 


Vf 




P2 


Ps 


Li 


1.0000 


1.0000 


1.0000 


1.0000 


-24.7916 


0.0000 


Al 


3.0000 


3.0000 


3.0000 


8.0000 


-23.1826 


0.0076 


Si 


4.0000 


4.0000 


4.0000 


4.0000 


-4.1684 


21.7115 





2.0000 


6.0000 


4.0000 


4.0000 


-3.5500 


37.5000 




Pv3 


Pv5 


Plp2 


Pboc'i 


Pboci 


Pbocb 


Li 


2.2989 


2.8103 


0.0000 


6.9107 


5.4409 


0.1973 


Al 


1.5000 


2.5791 


0.0000 


0.2500 


20.0000 


0.0000 


Si 


2.0754 


2.5791 


0.0000 


23.8188 


9.0751 


0.8381 


o 


2.9000 


2.9225 


0.4056 


0.7640 


3.5027 


0.0021 



TABLE A. VI: Bond order parameters. 



Bond 


Pbo,l 


Pbo,2 


Pbo,3 


PboA 


Pbo,5 


Pbo,6 


O- 


O 


5.5000 


1. 


.0000 


9.0000 


1 


.0000 


-0.1000 


0.6051 


Si- 


-0 


8.4790 


6. 


.0658 


28.8153 


1 


.0000 


-0.3000 


0.2131 


Si- 


-Si 


8.7229 


0. 


.0000 


7.1248 


1 


.0000 


-0.3000 


0.0392 


Al- 


-0 


6.1462 


0. 


.0000 


25.0000 


1 


.0000 


-0.3000 


0.1925 


Al- 


-Si 


10.0000 


0. 


.0000 


12.0000 


1 


.0000 


0.3000 


1.0000 


Al- 


-Al 


6.1608 


0. 


.0000 


14.3085 


1 


.0000 


-0.3000 


0.5154 


Li- 


-o 


7.8656 


0. 


.0000 


11.9965 


1 


.0000 


0.3000 


0.3228 


Li- 


-Si 


10.0000 


0. 


.0000 


12.0000 


1 


.0000 


0.3000 


1.0000 


Al- 


-Li 


10.0000 


0. 


.0000 


12.0000 


1 


.0000 


0.3000 


1.0000 


Li- 


-Li 


4.0000 


0. 


.0000 


12.0000 


1 


.0000 


0.3000 


0.6003 



TABLE A.III: Covalent radii [r%, rff, rjT in A] and Coulomb 
interaction parameters [rj (eV), x (eV) and 7 (A)]. 











Coulomb parameters 


Atom 


rZ 


r% 


rZ n 


V X 7 


Li 


1.6908 


-0.1000 


-1.0000 


11.0234 -3.2182 1.0000 


Al 


2.1967 


-1.6836 


-1.0000 


6.5000 -0.3343 0.4961 


Si 


2.1932 


1.2962 


-1.0000 


5.5558 4.2033 0.5947 


O 


1.2450 


1.0548 


0.9049 


8.3122 8.5000 1.0898 



TABLE A. VII: Off-diagonal bond parameters [Aj 
(kcal/mol), a (unitless)] and bond radii [R v dw, »*o s r o\ 
and rZ n (A)]. 



Bond 




RvdW 


a 


a 

ro 


rZ 


7T7T 


Si- 


O 


0.1836 


1.9157 


10.9070 


1.7073 


1.2375 


-1.0000 


Al- 


-O 


0.2017 


1.8458 


11.0700 


1.6009 


-1.0000 


-1.0000 


Al- 


-Si 


0.1000 


1.8500 


10.3237 


-1.0000 


-1.0000 


-1.0000 


Li- 


-0 


0.0790 


2.2000 


9.0491 


1.8165 


-1.0000 


1.0000 


Li- 


-Si 


0.0200 


1.5000 


10.0529 


-1.0000 


1.0000 


1.0000 


Li- 


-Al 


0.1146 


2.2000 


9.7537 


-1.0000 


1.0000 


1.0000 
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TABLE A. VIII: Valence angle parameters 









e , 

(deg.) 


Pvl 

(kcal/mol) 


Pv2 


PvA 


Pv7 


o- 


0- 





80.7324 


30.4554 


0.9953 


1.0783 


1.6310 


Si- 


-Si- 


-Si 


78.5339 


36.4328 


1.0067 


1.6608 


0.1694 


0- 


•Si- 


-Si 


86.3294 


18.3879 


5.8529 


1.2310 


1.7361 


0- 


•Si- 


-0 


79.5581 


34.9140 


1.0801 


2.2206 


0.1632 


Si- 


0- 


-Si 


82.3364 


4.7350 


1.3544 


1.0400 


1.4627 


0- 


O- 


-Si 


92.1207 


24.3937 


0.5000 


3.0000 


1.7208 


0- 


o- 


Al 


34.4326 


25.9544 


5.1239 


1.7141 


2.7500 


Al- 


-0- 


-Al 


20.7204 


13.4875 


4.0000 


1.4098 


0.6619 


0- 


•Al- 


-O 


59.5433 


20.0000 


4.0000 


2.0988 


3.0000 


0- 


•Li- 


-0 


60.0000 


0.0000 


1.0000 


1.0000 


1.0000 


0- 


O- 


Li 


81.6233 


30.0000 


2.0000 


1.0000 


1.0000 


Li- 


-0- 


-Li 


67.5247 


6.4512 


4.0000 


2.8079 


1.0000 


Al- 


-0- 


-Li 


50.9423 


7.0901 


3.9271 


2.5544 


1.0000 


Si- 


0- 


-Al 


18.0953 


5.3220 


4.0000 


1.0139 


1.0000 


Si- 


•o- 


-Li 


62.6634 


8.4441 


2.5120 


1.0000 


1.0000 



Appendix B: Calculation of elastic constants 

The elements of the elastic stiffness tensor Cijki for a 
and j3 eucryptite were computed within the framework 
of ReaxFF by calculating the second derivatives of strain 
energy density with respect to the strain components 89 



ijkl 



d 2 {E/V) 

dCijEkl 



(Bl) 



where E is the elastic energy stored in a domain of vol- 
ume V of the crystal subjected to homogeneous defor- 
mations. A similar approach has been employed earlier 
for computing the elastic constants of /3-eucryptite using 
DFT calculations (See Ref. 76). For sufficiently small 
strains, the total energy E of a crystal subjected to a 
general strain can be expressed as a Taylor series expan- 
sion truncated at the second order— 



E{V,e) 




El 

h3 



Vi Vj 



(B2) 

where the subscripts are cast in the Voigt notation (11=1, 
22=2, 33=3, 23=4, 31=5, and 12=6), % = 1 if i = 
1, 2, or 3 and r]i = 2 if i = 4, 5, or 6, E is the energy of 
the crystal volume Vq at equilibrium, cr^ are the elements 
of the stress tensor, and Sij is the Kronecker symbol. For 
the strains listed in Tables |B~T1 and iBTfl Eq. (|B2|) reduces 
to 



E(V,8) = E + V a (A 1 S + A 2 S 2 ), 



(B3) 



where A\ is related to stress components , and Ai is a 
linear combination of the elastic constants expressed 
in the Voigt notation. 

/3-eucryptite has five independent elastic constants 
namely, Cn, C12, C13, C33 and C44 due to the hexag- 
onal symmetry associated with its structured Table IB. II 



TABLE B.I: The strains used to calculate the five indepen- 
dent elastic constants of hexagonal /3-eucryptite (also used in 
Refs. HO and[7i). The relationship between A 2 in Eq. (|B3)| 
and dj are also provided. 



Strain parameters 
(unlisted ti — 0) 



Second-order coefficient 
A 2 in Eq. (|B3l> 



ei = £2 = 8 C11 + C12 

€1 — —62 = 5 Cll — C12 

ei = e 2 = ea = 8 Cii + C12 + 2C13 + C33/2 

£3=5 C33/2 

£5 = 5 2C44 



TABLE B.II: The strains used to calculate the seven inde- 
pendent elastic constants of rhombohedral a-eucryptite. The 
relationship between A2 in Eq. (|B3|I and dj are also provided. 



Strain parameters 
(unlisted ei = 0) 



Second-order coefficient 
A 2 in Eq. fB3) 



£1=5 
£3 = 5 
£1 = £2 = £3 = 8 
£l = -£2 =£3=5 
£l = £4 = 5 

£1 = £5 = 5 
£4 = 5 



C*n/2 

C33/2 

Cii + C12 + 2C13 + C33/2 
C11 — C12 + C33/2 
Cn/2 + 2C14 + 2C44 
Cn/2 + 2Ci5 + 2C44 
2C44 



lists the five different strains that we utilised to compute 
the elastic constants of /3-eucryptite along with the rela- 
tionship between the second-order coefficient A 2 and the 
elastic constants CV,- for each type of strain. 

On the other hand, a-eucryptite has a rhombohedral 
crystal structure and thereby, has seven independent elas- 
tic constants namely, Cn, C12, C13, C14, C15, C33 and 
C44^ The different strains used to compute these seven 
elastic constants and the relationships between A2 and 
Cij for each type of strain have been summarized in Ta- 
ble EHi 

For a given crystal, the total energy was computed for 
different values of 8 ranging from -2% to 2% using the 
LAMMPS2! implementation of ReaxFF. The calculated 
data were then fit to Eq. (|B3[) to extract the second- 
order coefficients A2 which were then used to evaluate 
the elastic constants through the relationships given in 
Tables ED and EUl 
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